# datos rasterrizados y espaciales
library(ncdf4)
library(raster)
library(sf)
library(CCAMLRGIS)
library(tmaptools)
# manipulacion
library(readr)
library(tidyverse)
library(patchwork)
library(data.table)
# analisis estadisticos post
library(easystats)
library(see)
# visualizacion
library(ggridges)
library(ncdf4)
# optimización
library(parallel)
#library(kableExtra)
library(devtools)
# datos precargados de hielo marino
#devtools::install_github("AustralianAntarcticDivision/raadtools")
library(raadtools)
#roots
library(here)
Este documento tiene como objetivo, en primera instancia, entender los datos de diversas variables ambientales que influyen en la dinámica poblacional del krill. En este caso se trabaja la lectura, manipulación y extracción de series temporales en toda la Antártica de las variables de interés, para despues centrar los datos en la zona de estudio, es decir, la Península Antártica. Una vez obtenidas las series temporales, se generan bases de datos que servirán para realizar test de correlación estacio temporal con variabes biologica-pesquera del krill, en este caso, tallas medias y tambien tallas de reclitamiento, todos estos indicadores siendo dependientes de la pesquería.
una vez entendidas estas relaciones, se ultilizaran las conculisiones y datos para incorporar en un modelo de evaluación integrada del krill en la PA:
.ncContributor_name: Walter N. Meier, Florence Fetterer, Ann Windnagel, J. Scott Stewart, Trey Stafford, Matt Fisher
This project was supported in part by a grant from the NOAA Climate Data Record Program. Production of original NASA Team and Bootstrap algorithm estimates supported by the NASA Polar Distributed Active Archive Center. The sea ice concentration algorithms were developed by Donald J. Cavalieri, Josefino C. Comiso, Claire L. Parkinson, and others at the NASA Goddard Space Flight Center in Greenbelt, MD. platform: NIMBUS-7 Sensor: SMMR > Scanning Multichannel Microwave Radiometer
Un mapa de que ilustra los cambios en la concentracion de hielo a través de los años.
spplot(ich,
colorkey = list(space = "bottom"),
scales = list(draw = TRUE),
main = "Índice Concentracion de Hielo",
fill='white',
names.attr=c("1978", "1979","1980", "1981", "1982",
"1983","1984", "1985", "1986",
"1987","1989", "1990", "1991",
"1992", "1993", "1994",
"1995","1996", "1997", "1998",
"1999", "2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019", "2020", "2021"))
Ahora preparo los datos como data frame.
extension <- extent(ich)
proyection <- crs(ich)
# ahora extraigo los valores por cada punto
stack_zona_ich <- raster::as.data.frame(ich, xy = TRUE)
stack_zona_ich <- na.omit(stack_zona_ich)
Ahora debo manipuñar la data para dejar formato largo table y luego extraer la serie de tiempo.
# Primero cambio nombre de las columnas
colnames(stack_zona_ich) <- c("lon", "lat", "1978",
"1979","1980", "1981", "1982",
"1983","1984", "1985", "1986",
"1987","1989", "1990", "1991",
"1992", "1993", "1994",
"1995","1996", "1997", "1998",
"1999", "2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019", "2020", "2021")
# ahora pivoteo los datos a lo largo
datsic <- stack_zona_ich %>%
pivot_longer(cols = c("1978", "1979","1980", "1981", "1982",
"1983","1984", "1985", "1986",
"1987","1989", "1990", "1991",
"1992", "1993", "1994",
"1995","1996", "1997", "1998",
"1999", "2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019", "2020", "2021",),
names_to = "ANO", values_to = "SEAICE")
Una vez ordenada los datos como un data.frame, puedo
obtener los registros agregados como serie de tiempo.
Un histograma de para ver la distribución de la variable y otros estadisticos adhoc
dsic <- datsic %>%
group_by(ANO) %>%
summarize(measic=mean(SEAICE))
# trend
we <- ggplot(dsic,
aes(x = as.double(ANO),
y = measic)) +
geom_line(col="blue")+
geom_point()+
stat_smooth(method = "loess", col="blue")+
theme_bw()+
theme( panel.grid.major = element_blank())+
xlim(1979, 2021)+
ylim(0.60, 0.75)+
ylab("Indice Concentración de Hielo")+
xlab("")
# histograma
weh <- ggplot(dsic,
aes(x= measic)) +
geom_histogram(col=2, fill=2,
bins = 50)+
theme_bw()
weh|we
.ncsetwd("~/DOCAS/Data/KrillLows/Data_All")
nc <- nc_open("antarctic_omi_si_extent_19930115_P20220328.nc")
# Obtener una lista de las variables en el archivo
attributes(nc$var)
## $names
## [1] "siextents_cglo" "siextents_oras" "siextents_foam" "siextents_glor"
## [5] "siextents_mean" "siextents_std"
attributes(nc$dim)
## $names
## [1] "time"
# Obtener una variable específica
var <- ncvar_get(nc, "siextents_glor")
time_var <- ncvar_get(nc, "time")
# Crear dtaframe
df <- data.table(var)
time_var_r <- as.POSIXct(time_var, origin = "1950-01-01 00:00:00", tz = "UTC")
df <- data.frame(time_var_r,var)
data <- df %>%
group_by(time_var_r) %>%
summarize(mean_var=mean(var))
#trend
we1<- ggplot(data , aes(x = time_var_r, y = mean_var)) +
geom_line()+
stat_smooth(method = "loess")+
theme_bw()
# histograma
weh1 <- ggplot(data,
aes(x= mean_var)) +
geom_histogram(col=4, fill=4,
bins = 50)+
theme_bw()
weh1|we1
Sin embargo estos datos solo corresponden a un año de registros, el
2022. Solo leí el archivo para ejemplificar la lectura de este tipo de
archivos.
The Sea Ice Index provides a quick look at Arctic- and Antarctic-wide changes in sea ice. It is a source for consistent, up-to-date sea ice extent and concentration images, in PNG format, and data values, in GeoTIFF and ASCII text files, from November 1978 to the present. Sea Ice Index images also depict trends and anomalies in ice cover calculated using a 30-year reference period of 1981 through 2010. The images and data are produced in a consistent way that makes the Index time-series appropriate for use when looking at long-term trends in sea ice cover. Both monthly and daily products are available. However, monthly products are better to use for long-term trend analysis because errors in the daily product tend to be averaged out in the monthly product and because day-to-day variations are often the result of short-term weather.
# Crear un vector con la lista de nombres de archivo
file_list <- list.files(path = "~/DOCAS/Data/KrillLows/01_Jan",
pattern = "*.tif", full.names = TRUE)
# Leer cada archivo en la lista con lappl
raster_lis <- lapply(file_list, raster)
# Unir los archivos en un único objeto "raster stack"
raster_stack <- stack(raster_lis)
extension <- extent(raster_stack)
proyeccion <- crs(raster_stack)
spplot(raster_stack,
colorkey = list(space = "bottom"),
scales = list(draw = TRUE),
main = "Concentracion de Hielo",
fill='white',
names.attr=c("1979","1980", "1981", "1982",
"1983","1984", "1985", "1986",
"1987","1989", "1990", "1991",
"1992", "1993", "1994",
"1995","1996", "1997", "1998",
"1999", "2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019", "2020", "2021", "2022"))
# Seleccionar la primera capa del "raster stack"
#raster_interes <- raster_stack[[1]]
# ahora extraigo los valores por cada punto
stack_zona_df <- raster::as.data.frame(raster_stack, xy = TRUE)
stack_zona_df <- na.omit(stack_zona_df)
#str(stack_zona_df)
Plot simples de la variable desde archivos .tif de un
año ejemplo (2001)
p2001 <- ggplot() +
geom_histogram(data = stack_zona_df, aes(S_201701_concentration_v3.0))+
theme_bw()+
xlab("")
p2001a <- ggplot() +
geom_raster(data = stack_zona_df, aes(x = x, y = y, fill = S_200101_concentration_v3.0))+
scale_fill_fermenter( palette = "Blues",
name="Sea Ice Concentracion 2001")+
theme_bw()+
xlab("")
p2001|p2001a
Ahora debo manipuñar la data para dejar formato largo table y luego extraer la serie de tiempo.
# Primero cambio nombre de las columnas
colnames(stack_zona_df) <- c("lon", "lat",
"1979","1980", "1981", "1982",
"1983","1984", "1985", "1986",
"1987","1989", "1990", "1991",
"1992", "1993", "1994",
"1995","1996", "1997", "1998",
"1999", "2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019", "2020", "2021", "2022")
# ahora pivoteo los datos a lo largo
datgeo2 <- stack_zona_df %>%
pivot_longer(cols = c("1979","1980", "1981", "1982",
"1983","1984", "1985", "1986",
"1987","1989", "1990", "1991",
"1992", "1993", "1994",
"1995","1996", "1997", "1998",
"1999", "2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019", "2020", "2021", "2022"),
names_to = "ANO", values_to = "SEAICE")
# preparo la data para transformar a lat long
dat_utm <- datgeo2[, c("lon", "lat")]
# Crea un objeto de proyección para UTM Zone 30N
utm30 <- CRS("+proj=utm +zone=30 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
# Transforma las coordenadas UTM a latitud y longitud
df_latlong <- spTransform(SpatialPoints(dat_utm, proj4string=utm30), CRS("+proj=longlat +datum=WGS84"))
# Extrae las coordenadas latitud y longitud del objeto transformado
df_latlong <- as.data.frame(df_latlong@coords)
# Renombra las columnas de las coordenadas
colnames(df_latlong) <- c("long", "lat")
datgeo4 <- cbind(df_latlong, datgeo2[,c(3,4)])
Compruebo la tendencia general de los datos a través de la media de concentracion de hielo en toda la zona
data=datgeo4 %>%
group_by(ANO) %>%
summarize(mean_var=mean(SEAICE))
# un grafico simple
as <- ggplot(data,aes(y=mean_var, x=as.double(ANO))) +
geom_point()+
stat_smooth(method = "loess")+
theme_light()
as
1 .Clhorophila (1997-2022) Global Ocean Colour (Copernicus-GlobColour), Bio-Geo-Chemical, L4 (monthly and interpolated) from Satellite Observations Overview
The biogeochemical hindcast for global ocean is produced at Mercator-Ocean (Toulouse. France). It provides 3D biogeochemical fields for the time period 1993-2019 at 1/4 degree and on 75 vertical levels. It uses PISCES biogeochemical model (available on the NEMO modelling platform). No data assimilation in this product.
Latest NEMO version (v3.6_STABLE) Forcings: GLORYS2V4-FREE ocean physics produced at Mercator-Ocean and ERA-Interim atmosphere produced at ECMWF at a daily frequency Outputs: Daily (chlorophyll. nitrate. phosphate. silicate. dissolved oxygen. primary production) and monthly (chlorophyll. nitrate. phosphate. silicate. dissolved oxygen. primary production. iron. phytoplankton in carbon) 3D mean fields interpolated on a standard regular grid in NetCDF format. The simulation is performed once and for all. Initial conditions: World Ocean Atlas 2013 for nitrate. phosphate. silicate and dissolved oxygen. GLODAPv2 for DIC and Alkalinity. and climatological model outputs for Iron and DOC Quality/Accuracy/Calibration information: See the related Quality Information Document DOI (product): https://doi.org/10.48670/moi-00019
Full name: Global ocean biogeochemistry hindcast Product ID : GLOBAL_MULTIYEAR_BGC_001_029 DOI:10.48670/moi-00019 Source :Numerical models Spatial extent : Global OceanLat -90° to 90°Lon -180° to 180° Spatial resolution: 0.25° × 0.25° Temporal extent: 1 Jan 1993 to 31 Dec 2020 Temporal resolution: Daily-Monthly Elevation (depth) levels: 75 Processing level: Level 4 Variables: Cell thickness Mass concentration of chlorophyll a in sea water (CHL) Model level number at sea floorMole concentration of dissolved iron in sea water (FE) Mole concentration of dissolved molecular oxygen in sea water (O2) Mole concentration of nitrate in sea water (NO3) Mole concentration of phosphate in sea water (PO4) Mole concentration of phytoplankton expressed as carbon in sea water (PHYC) Mole concentration of silicate in sea water (SI) Net primary production of biomass expressed as carbon per unit volume in sea water (PP) Sea floor depth below geoidSea water ph reported on total scale (pH) Surface partial pressure of carbon dioxide in sea water (spCO2) Feature type: Grid Blue markets: Conservation & biodiversityClimate & adaptationPolicy & governanceScience & innovationMarine food Projection: ETRS89 (EPSG 4258) Data assimilation: None Update frequency: Weekly Format: NetCDF-4 Originating centre: Mercator Océan International filename: cmems_mod.nc
Primero inspeccionamos el archivo y sus atributos.
chl <- nc_open(here("Data_TSM_Chla","cmems_mod.nc"))
# Obtener una lista de las variables en el archivo
attributes(chl$var)
## $names
## [1] "chl"
attributes(chl$dim)
## $names
## [1] "time" "depth" "latitude" "longitude"
En este caso, solo tenemos una variable dasociada al raster llamada
chl
Ahora leo los datos por fecha para manipular mejor el archivo
.nc con la función brick(). (A RasterBrick is
a multi-layer raster object. They are typically created from a
multi-layer (band) file; but they can also exist entirely in memory.
They are similar to a RasterStack (that can be created with stack), but
processing time should be shorter when using a RasterBrick. Yet they are
less flexible as they can only point to a single file.)
# Abrir el archivo .nc
setwd("~/DOCAS/Data/KrillLows/Data_All")
raster_nc <- brick("cmems_mod.nc", varname = "chl")
# Obtener la variable temporal
dates <- getZ(raster_nc)
years <- as.integer(format(as.Date(dates), "%Y"))
month <- as.integer(format(as.Date(dates), "%m"))
unique(month)
## [1] 12 1 2 3 4 5 6 7 8 9 10 11
chla <- stackApply(raster_nc, years, fun=mean)
extension <- extent(chla)
proyeccion <- crs(chla)
# ahora extraigo los valores por cada punto
stack_zona_chl <- raster::as.data.frame(chla, xy = TRUE)
stack_zona_chl <- na.omit(stack_zona_chl)
#str(stack_zona_chl)
Saco un plot del comportamiento de la variable en todos los años
spplot(chla,
colorkey = list(space = "bottom"),
scales = list(draw = TRUE),
main = "Chorophyla a",
fill='white',
names.attr=c("2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019", "2020"))
Miro la distribución de la variable en el año 2000
ggplot() +
geom_histogram(data = stack_zona_chl,
aes(index_2020),
bins= 100)+
theme_bw()
Ahora ploteo los datos del 2000
ggplot() +
geom_raster(data = stack_zona_chl,
aes(x = x, y = y, fill = index_2000))+
theme_bw()+
scale_fill_fermenter(palette = "BrBG")
Ahora debo manipuñar la data para dejar formato largo table y luego
extraer la serie de tiempo.
# Primero cambio nombre de las columnas
colnames(stack_zona_chl) <- c("lon", "lat",
"2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019", "2020")
# ahora pivoteo los datos a lo largo
datchl <- stack_zona_chl %>%
pivot_longer(cols = c( "2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019", "2020"),
names_to = "ANO", values_to = "CHLa")
Compruebo la tendencia general de los datos a traves de la media de concentracion de hielo
daa=datchl %>%
group_by(ANO) %>%
summarize(mean_var=mean(CHLa))
# un grafico simple
aschl <- ggplot(daa,aes(y=mean_var, x=as.double(ANO))) +
geom_point(color="2")+
stat_smooth(method = "lm", col=2)+
theme_light()+
xlab("")
aschl
ERA5 is the fifth generation ECMWF reanalysis for the global climate and weather for the past 4 to 7 decades. Currently data is available from 1950, with Climate Data Store entries for 1950-1978 (preliminary back extension) and from 1959 onwards (final release plus timely updates, this page). ERA5 replaces the ERA-Interim reanalysis.
Reanalysis combines model data with observations from across the world into a globally complete and consistent dataset using the laws of physics. This principle, called data assimilation, is based on the method used by numerical weather prediction centres, where every so many hours (12 hours at ECMWF) a previous forecast is combined with newly available observations in an optimal way to produce a new best estimate of the state of the atmosphere, called analysis, from which an updated, improved forecast is issued. Reanalysis works in the same way, but at reduced resolution to allow for the provision of a dataset spanning back several decades. Reanalysis does not have the constraint of issuing timely forecasts, so there is more time to collect observations, and when going further back in time, to allow for the ingestion of improved versions of the original observations, which all benefit the quality of the reanalysis product.
ERA5 provides hourly estimates for a large number of atmospheric, ocean-wave and land-surface quantities. An uncertainty estimate is sampled by an underlying 10-member ensemble at three-hourly intervals. Ensemble mean and spread have been pre-computed for convenience. Such uncertainty estimates are closely related to the information content of the available observing system which has evolved considerably over time. They also indicate flow-dependent sensitive areas. To facilitate many climate applications, monthly-mean averages have been pre-calculated too, though monthly means are not available for the ensemble mean and spread.
ERA5 is updated daily with a latency of about 5 days (monthly means are available around the 6th of each month). In case that serious flaws are detected in this early release (called ERA5T), this data could be different from the final release 2 to 3 months later. In case that this occurs users are notified.
The data set presented here is a regridded subset of the full ERA5 data set on native resolution. It is online on spinning disk, which should ensure fast and easy access. It should satisfy the requirements for most common applications.
An overview of all ERA5 datasets can be found in this article. Information on access to ERA5 data on native resolution is provided in these guidelines.
Data has been regridded to a regular lat-lon grid of 0.25 degrees for the reanalysis and 0.5 degrees for the uncertainty estimate (0.5 and 1 degree respectively for ocean waves). There are four main sub sets: hourly and monthly products, both on pressure levels (upper air fields) and single levels (atmospheric, ocean-wave and land surface quantities).
This file have “ERA5 monthly mean data on single levels from 1991 to present”.
tsm <- nc_open("~/DOCAS/Data/KrillLows/Data_TSM_Chla/tsm_1995_2022_PA.nc")
# Obtener una lista de las variables en el archivo
attributes(tsm$var)
## $names
## [1] "sst"
attributes(tsm$dim)
## $names
## [1] "longitude" "latitude" "expver" "time"
Ahora leo los datos por fecha para manipular mejor
# Abrir el archivo .nc
setwd("~/DOCAS/Data/KrillLows/Data_All")
raster_ts <- brick("tsm_1995_2022_PA.nc", varname = "sst")
# Obtener la variable temporal
dates <- getZ(raster_ts)
years <- as.integer(format(as.Date(dates), "%Y"))
#month <- as.integer(format(as.Date(dates), "%m"))
unique(years)
## [1] 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005
## [16] 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020
## [31] 2021 2022
#unique(month)
sstm <- stackApply(raster_ts, years, fun=mean)
extension <- extent(sstm)
proyeccion <- crs(sstm)
# ahora extraigo los valores por cada punto
stack_zona_tsm <- raster::as.data.frame(sstm, xy = TRUE)
stack_zona_tsm <- na.omit(stack_zona_tsm)
str(stack_zona_tsm)
## 'data.frame': 32012 obs. of 34 variables:
## $ x : num -100 -99.8 -99.5 -99.2 -99 ...
## $ y : num -50 -50 -50 -50 -50 -50 -50 -50 -50 -50 ...
## $ index_1991: num 281 281 281 281 281 ...
## $ index_1992: num 281 281 281 281 281 ...
## $ index_1993: num 281 281 281 281 281 ...
## $ index_1994: num 281 281 281 281 281 ...
## $ index_1995: num 281 281 281 281 281 ...
## $ index_1996: num 281 281 281 281 281 ...
## $ index_1997: num 281 281 281 281 281 ...
## $ index_1998: num 282 282 282 282 282 ...
## $ index_1999: num 281 281 281 281 281 ...
## $ index_2000: num 281 281 281 281 281 ...
## $ index_2001: num 281 281 281 281 281 ...
## $ index_2002: num 281 281 281 281 281 ...
## $ index_2003: num 281 281 281 281 281 ...
## $ index_2004: num 281 281 281 281 281 ...
## $ index_2005: num 281 281 281 281 281 ...
## $ index_2006: num 281 281 281 281 281 ...
## $ index_2007: num 281 281 281 281 281 ...
## $ index_2008: num 281 281 281 281 281 ...
## $ index_2009: num 281 280 280 280 280 ...
## $ index_2010: num 281 281 281 281 281 ...
## $ index_2011: num 280 280 280 280 280 ...
## $ index_2012: num 281 281 281 281 281 ...
## $ index_2013: num 281 281 281 281 281 ...
## $ index_2014: num 281 281 281 281 281 ...
## $ index_2015: num 281 281 281 281 281 ...
## $ index_2016: num 282 282 282 282 282 ...
## $ index_2017: num 281 281 281 281 281 ...
## $ index_2018: num 281 281 281 281 280 ...
## $ index_2019: num 281 281 281 281 281 ...
## $ index_2020: num 281 281 281 281 281 ...
## $ index_2021: num 281 281 281 281 281 ...
## $ index_2022: num 281 281 281 281 281 ...
## - attr(*, "na.action")= 'omit' Named int [1:22889] 102 103 104 105 106 107 108 109 110 111 ...
## ..- attr(*, "names")= chr [1:22889] "102" "103" "104" "105" ...
Un plot de los datos a traves de los años
spplot(sstm,
colorkey = list(space = "bottom"),
scales = list(draw = TRUE),
main = "TSM",
fill='white')
Ahora ploteo los datos del 2000
ggplot() +
geom_raster(data = stack_zona_tsm,
aes(x = lon, y = lat, fill = stack_zona_tsm$`2000`))+
theme_bw()+
scale_fill_fermenter( palette = "Oranges",
n.breaks = 7,
name="TSM 2000")+
ylim(-75,-51)+
xlim(-100, -30)
plot(sstm, col = colorRampPalette(c("blue", "green", "yellow", "red"))(10), main = "Mapa de la variable de interés", ylim=c(-70,-50))
Ahora debo manipuñar la data para dejar formato largo table y luego
extraer la serie de tiempo.
# Primero cambio nombre de las columnas
colnames(stack_zona_tsm) <- c("lon", "lat",
"1991","1992", "1993", "1994",
"1995","1996", "1997", "1998",
"1999", "2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019", "2020",
"2021", "2022")
# ahora pivoteo los datos a lo largo
dattsm <- stack_zona_tsm %>%
pivot_longer(cols = c("1991","1992", "1993", "1994",
"1995","1996", "1997", "1998",
"1999", "2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019", "2020",
"2021", "2022"),
names_to = "ANO", values_to = "TSM")
Compruebo la tendencia general de los datos a traves de la media TSM
Distribución de la variable en frecuencia y tiempo
tsmh <- ggplot() +
geom_histogram(data = dattsm,
aes(TSM),
bins= 100)+
theme_bw()
datsm1=dattsm %>%
group_by(ANO) %>%
summarize(mean_var=mean(TSM))
# un grafico simple
aschl <- ggplot(datsm1,aes(y=mean_var,
x=as.double(ANO))) +
geom_point(color="2")+
stat_smooth(method = "loess", col=2)+
theme_light()+
xlab("")
tsmh |aschl
# AGRUPACION DE VARIABLES EN EL ESPACIO
En este paso identificamos los distintos niveles de agregacion
espacial en klos cuales se realizará el analisis de las tendencias de
las variabless ambientales (Chla, Sea Ice cover y TSM). En este caso
usaremos grillas de 1x1 grados, las áreas de manejo espacial de CCAMLR y
el bolque 48.1. Primero identifico todos los objetos
RasterStack definidos previamente.
SIC = ich Datos de Hielo de archivos .nc
SIC =raster_stack Datos dde hielo de archivos
.tif (1979-2022) Chla = chla datos de
clorofila en archivo .nc TSM = sstm Datos de
TSM en archivo unico .ncdesde el año 1991 al 2022
Identificar tendencias temporales por cada SSSMU para tener un mejor detalle del comportamiento de la variable a traves del tiempo en el espacio definido.
Defino las subareas a trabajar con la librería
CCAMLRGIS
# con SSMU
ssmu <- load_SSMUs()
ssmu481 <- subset(ssmu[(2:7),])
ssmu481a <- st_as_sf(ssmu481)
ssmu481aa = st_transform(ssmu481a, "+proj=latlong +ellps=WGS84")
# con Statistical Areas
suba <- load_ASDs()
suba1 <- subset(suba[(3),])
suba1a<- st_as_sf(suba1)
suba1aa = st_transform(suba1a, "+proj=latlong +ellps=WGS84")
# Cargo linea de costa
coast <- load_Coastline()
coast1<- st_as_sf(coast)
coast2 = st_transform(coast1, "+proj=latlong +ellps=WGS84")
# y testeo el mapa
ssmap <- ggplot(ssmu481aa)+
geom_sf(data = coast2, colour="black", fill=NA)+
geom_sf(data= suba1aa, fill=NA)+
geom_sf(aes(fill=ssmu481aa$GAR_Short_Label))+
scale_fill_viridis_d(option = "E")+
geom_sf_label(aes(label = ssmu481aa$GAR_Short_Label))+
labs(fill = "SSMU")+
ylim(-70, -60)+
xlim(-70, -50)+
theme_bw()
ssmap
De forma experimental hago una grilla sobre las SSMU.
Grid<- ssmu481aa %>% #pm481 es el plot base original linea 481
sf::st_make_grid(cellsize = c(0.5,0.25)) %>% # para que quede cuadrada
sf::st_cast("MULTIPOLYGON") %>%
sf::st_sf() %>% # objeto en spatial feature
dplyr::mutate(cellid = row_number())
# Corto la grilla dentro de las SSMU
gridcrop <- crop_shape(Grid, ssmu481aa, polygon = TRUE)
# Pruebo el mapa
ggplot() +
geom_sf(data = gridcrop,
color="grey",
fill=NA,
size=0.3) +
geom_sf(data = ssmu481aa, fill = NA)+
geom_sf(data= suba1aa, fill=NA)+
geom_sf(data = coast2, colour="black", fill=NA)+
theme_bw()+
ylim(-70, -60)+
xlim(-70, -50)
Extraigo los datos previamente trabajados del objeto
raster_stack
extich<-raster::extract(raster_stack, ssmu481aa,
method="bilinear",
weights=FALSE,
buffer=NULL,
small=FALSE,
cellnumbers=FALSE,
fun=mean,
na.rm=TRUE,
df=TRUE,
factors=FALSE,
sp=TRUE)
Ahora lo ordeno como data.frame
coords<-as.data.frame(coordinates(extich))
daich<-data.frame(coords,extich)
# Primero cambio nombre de las columnas
daich1 <- daich %>%
select(1, 2, 5, 18:60)
colnames(daich1) <- c("LON", "LAT", "SSMU",
"1979","1980", "1981", "1982",
"1983","1984", "1985", "1986",
"1987","1989", "1990", "1991",
"1992", "1993", "1994",
"1995","1996", "1997", "1998",
"1999", "2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019",
"2020", "2021", "2022")
# ahora pivoteo los datos a lo largo
daich2 <- daich1 %>%
pivot_longer(cols = c("1979","1980", "1981", "1982",
"1983","1984", "1985", "1986",
"1987","1989", "1990", "1991",
"1992", "1993", "1994",
"1995","1996", "1997", "1998",
"1999", "2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019",
"2020", "2021", "2022"),
names_to = "ANO", values_to = "SIC")
Plot simples de la variable SIC en la zona de estudio
hissic <- ggplot() +
geom_histogram(data = daich2, aes(SIC, fill=SSMU))+
facet_wrap(~SSMU)+
scale_fill_bluebrown_d()+
theme_bw()
hissic
Grafico de la variable en el tiempo por cada SSMU
as <- ggplot(daich2,
aes(y=SIC, x=as.double(ANO))) +
geom_point()+
stat_smooth(method = "loess")+
facet_wrap(~SSMU, scales="free_y")+
theme_bw()+
xlab("")+
ylab("Sea Ice Concentration Mean")
as
Ahora un grafico por todas las SSMU
astot <- ggplot(daich2 %>%
group_by(ANO) %>%
summarize(meani=mean(SIC)),
aes(y=meani, x=as.double(ANO))) +
geom_point()+
stat_smooth(method = "loess")+
theme_bw()+
xlab("")+
ylab("Sea Ice Concentration Mean")
astot
### CHL a
Extraigo los datos previamente trabajados del objeto
chla
extchl<-raster::extract(chla, ssmu481aa,
method="bilinear",
weights=FALSE,
buffer=NULL,
small=FALSE,
cellnumbers=FALSE,
fun=mean,
na.rm=TRUE,
df=TRUE,
factors=FALSE,
sp=TRUE)
Ahora lo ordeno como data.frame
coordschl<-as.data.frame(coordinates(extchl))
dachl<-data.frame(coordschl,extchl)
# Primero cambio nombre de las columnas
dachl1 <- dachl %>%
select(1, 2, 5, 18:38)
colnames(dachl1) <- c("LON", "LAT", "SSMU",
"2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019",
"2020")
# ahora pivoteo los datos a lo largo
dachl2 <- dachl1 %>%
pivot_longer(cols = c("2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019",
"2020"),
names_to = "ANO", values_to = "CHL")
Plot del comportaniento de la variable por SSMU
hischl <- ggplot() +
geom_histogram(data = dachl2, aes(CHL, fill=SSMU))+
facet_wrap(~SSMU)+
scale_fill_bluebrown_d()+
theme_bw()
hischl
Tendencia temporal de la variable por SSMU
aschl <- ggplot(dachl2,
aes(y=CHL, x=as.double(ANO))) +
geom_point(col="#386641")+
geom_line(col="#386641")+
stat_smooth(method = "loess", col="#80ed99")+
facet_wrap(~SSMU, scales="free_y")+
theme_bw()+
xlab("")+
ylab("Chl a Concentration")
aschl
Ahora un grafico por todas las SSSMU
chtot <- ggplot(dachl2 %>%
group_by(ANO) %>%
summarize(meani=mean(CHL)),
aes(y=meani, x=as.double(ANO))) +
geom_point(col="#386641")+
stat_smooth(method = "loess", col="#80ed99")+
theme_bw()+
xlab("")+
ylab("Chl a Concentration")
chtot
Una estadistica del data frame
skimr::skim(dachl2)
| Name | dachl2 |
| Number of rows | 126 |
| Number of columns | 5 |
| _______________________ | |
| Column type frequency: | |
| character | 2 |
| numeric | 3 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| SSMU | 0 | 1 | 25 | 34 | 0 | 6 | 0 |
| ANO | 0 | 1 | 4 | 4 | 0 | 21 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| LON | 0 | 1 | -59.36 | 3.00 | -63.85 | -61.51 | -59.51 | -57.10 | -54.69 | ▃▇▃▃▃ |
| LAT | 0 | 1 | -62.59 | 1.05 | -64.37 | -63.33 | -62.48 | -61.69 | -61.19 | ▃▃▃▃▇ |
| CHL | 0 | 1 | 1.18 | 0.32 | 0.81 | 0.94 | 1.03 | 1.40 | 2.49 | ▇▃▂▁▁ |
Extraigo los datos previamente trabajados del objeto
sstm
extsm<-raster::extract(sstm, ssmu481aa,
method="bilinear",
weights=FALSE,
buffer=NULL,
small=FALSE,
cellnumbers=FALSE,
fun=mean,
na.rm=TRUE,
df=TRUE,
factors=FALSE,
sp=TRUE)
Ahora lo ordeno como data.frame
coords<-as.data.frame(coordinates(extsm))
datsm<-data.frame(coords,extsm)
# Primero cambio nombre de las columnas
datsm1 <- datsm %>%
select(1, 2, 5, 18:49)
colnames(datsm1) <- c("LON", "LAT", "SSMU",
"1991","1992", "1993", "1994",
"1995","1996", "1997", "1998",
"1999", "2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019",
"2020", "2021", "2022")
# ahora pivoteo los datos a lo largo
datsm2 <- datsm1 %>%
pivot_longer(cols = c("1991",
"1992", "1993", "1994",
"1995","1996", "1997", "1998",
"1999", "2000", "2001", "2002", "2003",
"2004", "2005", "2006", "2007", "2008", "2009",
"2010","2011", "2012", "2013", "2014", "2015",
"2016", "2017" , "2018", "2019",
"2020", "2021", "2022"),
names_to = "ANO", values_to = "TSM")
Plot simples de TSM
histsm <- ggplot() +
geom_histogram(data = datsm2, aes(TSM, fill=SSMU))+
facet_wrap(~SSMU)+
scale_fill_viridis_d(option="A")+
theme_bw()
histsm
Tendencia temporal del TSM por SSMU
atsm <- ggplot(datsm2,
aes(y=TSM, x=as.double(ANO))) +
geom_point(col="#ffb3c1")+
stat_smooth(method = "loess", col="#a4133c")+
facet_wrap(~SSMU, scales="free_y")+
theme_bw()+
xlab("")+
ylab("TSM Mean")
atsm
Ahora un grafico por todas las SSSMU
tstot <- ggplot(datsm2 %>%
group_by(ANO) %>%
summarize(meani=mean(TSM)),
aes(y=meani, x=as.double(ANO))) +
geom_point(col="#ffb3c1")+
stat_smooth(method = "loess", col="#a4133c")+
theme_bw()+
xlab("")+
ylab("TSM")
tstot
De forma experimental hago una grilla sobre la 48.1
gsic <- st_as_sf(datgeo5, coords = c("long", "lat"),
crs = "+proj=latlong +ellps=WGS84")
Grid2<- suba1aa %>% #pm481 es el plot base original linea 481
sf::st_make_grid(cellsize = c(0.6,0.3)) %>% # para que quede cuadrada
sf::st_cast("MULTIPOLYGON") %>%
sf::st_sf() %>% # objeto en spatial feature
dplyr::mutate(cellid = row_number())
# Corto la grilla dentro de las SSMU
gridcrop1 <- crop_shape(Grid2, suba1aa, polygon = TRUE)
#
resulsic <- Grid2 %>%
st_join(gsic) %>%
group_by(cellid)
Mapas
# Pruebo el mapa
machl <- ggplot() +
geom_sf(data=resultsic %>%
drop_na(), aes(fill = SEAICE),
alpha=0.7, color=NA) +
scale_fill_see_c(palette = "light")+
geom_sf(data = gridcrop1,
color=NA,
fill=NA) +
geom_sf(data = ssmu481aa, fill = NA)+
geom_sf(data= coast1, fill="#f9dcc4")+
geom_sf(data= suba1aa, fill=NA)+
xlab(expression(paste(Longitude^o,~'O'))) +
ylab(expression(paste(Latitude^o,~'S')))+
theme_bw()+
facet_wrap(~ANO)+
theme(panel.background = element_rect(fill = 'aliceblue'),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())+
ylim(-70, -60)+
xlim(-70, -50)
machl
De forma experimental hago una grilla sobre la 48.1
gtsm <- st_as_sf(dattsm , coords = c("lon", "lat"),
crs = "+proj=latlong +ellps=WGS84")
Grid2<- suba1aa %>% #pm481 es el plot base original linea 481
sf::st_make_grid(cellsize = c(0.6,0.3)) %>% # para que quede cuadrada
sf::st_cast("MULTIPOLYGON") %>%
sf::st_sf() %>% # objeto en spatial feature
dplyr::mutate(cellid = row_number())
# Corto la grilla dentro de las SSMU
gridcrop1 <- crop_shape(Grid2, suba1aa, polygon = TRUE)
#
result3 <- Grid2 %>%
st_join(gtsm) %>%
group_by(cellid)
Mapas
# Pruebo el mapa
matsm <- ggplot() +
geom_sf(data=result3 %>%
drop_na(), aes(fill = TSM),
alpha=0.7, color=NA) +
scale_fill_gradient(low = "#fff0f3",
high = "#d00000")+
geom_sf(data = gridcrop1,
color=NA,
fill=NA) +
geom_sf(data = ssmu481aa, fill = NA)+
geom_sf(data= coast1, fill="#f9dcc4")+
xlab(expression(paste(Longitude^o,~'O'))) +
ylab(expression(paste(Latitude^o,~'S')))+
theme_bw()+
facet_wrap(~ANO)+
theme(panel.background = element_rect(fill = 'aliceblue'),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())+
ylim(-70, -60)+
xlim(-70, -50)
matsm
De forma experimental hago una grilla sobre la 48.1
gchl <- st_as_sf(datchl, coords = c("lon", "lat"),
crs = "+proj=latlong +ellps=WGS84")
Grid2<- suba1aa %>% #pm481 es el plot base original linea 481
sf::st_make_grid(cellsize = c(0.6,0.3)) %>% # para que quede cuadrada
sf::st_cast("MULTIPOLYGON") %>%
sf::st_sf() %>% # objeto en spatial feature
dplyr::mutate(cellid = row_number())
# Corto la grilla dentro de las SSMU
gridcrop1 <- crop_shape(Grid2, suba1aa, polygon = TRUE)
#
resultchl <- Grid2 %>%
st_join(gchl) %>%
group_by(cellid)
Mapas
# Pruebo el mapa
machl <- ggplot() +
geom_sf(data=resultchl %>%
drop_na(), aes(fill = CHLa),
alpha=0.7, color=NA) +
scale_fill_distiller(palette = 5,
n.breaks = 5,
direction = -1)+
geom_sf(data = gridcrop1,
color=NA,
fill=NA) +
geom_sf(data = ssmu481aa, col="red", fill = NA)+
geom_sf(data= coast1, fill="#f9dcc4")+
geom_sf(data= suba1aa, fill=NA)+
xlab(expression(paste(Longitude^o,~'O'))) +
ylab(expression(paste(Latitude^o,~'S')))+
theme_bw()+
facet_wrap(~ANO)+
theme(panel.background = element_rect(fill = 'aliceblue'),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())+
ylim(-70, -60)+
xlim(-70, -50)
machl
Ahora se deben preparar las vbase de trabajo definitivas de las tres variables ambientales en cuestión; TSM, Concentración de Hielo y Chl-a
Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J., Nicolas, J., Peubey, C., Radu, R., Rozum, I., Schepers, D., Simmons, A., Soci, C., Dee, D., Thépaut, J-N. (2019): ERA5 monthly averaged data on single levels from 1959 to present. Copernicus Climate Change Service (C3S) Climate Data Store (CDS). (Accessed on < DD-MMM-YYYY >), 10.24381/cds.f17050d7